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1. Introduction 

Molecules with repeating subunits that spontaneously form helical shapes are ubiquitous 
in nature. A classical example is given by a-helices, which are an important motif in 
the secondary structure of proteins pQ. Other examples of naturally occurring helical 
structures include intermediate filament proteins keratin and vimentin p], the myosin 
and kinesin families of motor proteins [2], tubulin microtubules [3], RecA and RadSl 
filaments on DNA [1], flagella and pili [5], and others. Artificial polymer structures that 
spontaneously take helical forms have also been obtained, see for example [6l [7]. The 
spontaneous emergence and persistence of helical shapes has been the subject of several 
recent studies. For progress in the description of helical shapes based on elastic rod 
theory, see [H] • The spontaneous formation of helical structures under compaction of an 
elastic rod has also been shown, by using a Lennard- Jones-like potential whose repulsive 
core eliminates self- intersections [9]. 

Molecular shapes arise through an interplay between elastic forces caused by 
bending and twisting of the molecular bonds, and long- and short-range forces (such 
as electrostatic and van der Waals forces) between individual atoms comprising the 
molecule. For elastic interactions only, Kirchoff 's rod theory [101 HI] has been used 
to model molecular shapes with notable success [121 [ISl [El IISI [161 HI [HI- The 
recent work on non-local interactions has considered helical molecules as rods with 
charges concentrated along the axis of the rod only [HI [TH [H [19]. Other authors 
used energy optimization with short-range repulsion only, without initial assumptions 
on the molecular shapes, and yet helical shapes showed surprising persistence [9l [20] , 
either for the whole molecule, or in parts of it. The crucial question is whether such 
helical shapes persist for more complex interactions or geometries of the molecules. For 
example, more accurate models of molecules could include charges that are attached at 
a certain distance from the axis, as shown in Fig. [l| Interatomic forces are caused by 
some potential - or combination of potentials - that depends explicitly on the geometry 
of the molecule and the Euclidean distance d between two points s and s' along the 
curve representing the axis of the molecule. Equations of motion for such molecules 
were derived recently using nonlocal extensions of exact geometric rod theory [^ [22] • 
However, the question of existence of helical shapes for more complex interactions 
between the parts of the molecule and the methods for their computations has remained 
open. In this paper, we demonstrate a fast and efficient method for finding such states, 
and in addition show how we can achieve a complete classification of all helical states. We 
also show that the geometric approach allows exact computation of dispersion relations 
and linear instability growth rates. 

Plan of the paper. Section |2] derives the equations for nonlocally interacting rods of 
rigid charge configurations, called bouquets. These bouquets are fixed configurations of 
charges, rigidly attached to the molecule's central axis. The molecule is comprised 
of a long repeating chain of such bouquets. The equations are formulated on the 



Helical states of nonlocally interacting molecules 



3 



special Euclidean group, SE{?)), consisting of three-dimensional spatial rotations and 
translations, which naturally generate helical configurations. Section [3] demonstrates 
how to find helical stationary shapes by minimizing the energy for bouquets positioned 
along an infinitely long molecule with charges that are off-axis. As shown in Section 
|4| even this simple molecule, when deformed into helical configurations, shows quite 
a complex and intriguing energy landscape. Section |5] briefiy discusses how to extend 
the ideas presented above to include more general shapes. In particular, we consider 
a molecular shape we refer to as a 2-helix. Here, we define an n-helix as a molecule 
that is a true helix consisting of repeating groups of n bouquets along the axis. For 
a 2-helix, one optimizes the energy over both the shape of the helix and the relative 
configuration of two bouquets. One may treat n-helices by optimizing configurations 
of n bouquets using the same methods. After revealing the complexity of the energy 
landscape, even when limited to helical shapes, and the large number of energy minima 
in that "helical universe", it is natural to pose the question of linear stability and 
selection of the stationary shapes. By using an SE{?>) method that naturally accounts 
for the helical geometry, in Section [6] we investigate the linear stability of a general 
helical polymer that possesses torque-inducing non-local self-interactions. As far as 
we know, no such work has been undertaken previously. The geometric approach in 
Section |6] allows us to find the exact dispersion relation for the stability of the helical 
shapes with an arbitrary interaction potential. Of course, for particular applications 
it is important to explicitly compute the linearization operators. This is accomplished 
in Section [7] using geometric methods to compute derivatives of the potential energy. 
Section |8] computes the numerical stability for the particular example of a linear polymer 
comprising a charged rod with repeated configuration of two equal and opposite charges 
that interact through a screened electrostatic and Lennard- Jones potential. For the 
sake of simplicity, we concentrate on the stability of a polymer that is perfectly straight 
in its unstressed configuration. Such a polymer is neutrally stable in the absence of 
the nonlocal interactions. It is therefore interesting that nonlocal terms can induce 
instability that causes the molecule to deform into helical conformations. Physically, 
this instability is connected to the tendency of the rod to minimize its energy and 
properly align the dipole moments of each bouquet by twisting, as was already seen in 
the minimum energy calculation shown in Section |4| 

We shall note that these results are difficult to achieve in the traditional Kirchhoff- 
based approach, as the non-local interactions depend on the relative distance and 
orientation of the charges at different points on the rod. In these traditional methods, 
the equations of motion are written in a coordinate system that is moving with the 
rod, and changing with both position on the rod s and time t. Thus, one has to 
write equations of motion, i.e. calculate the momenta and forces, in a frame attached 
to the rod that is moving and rotating in a non-inertial fashion. The main difficulty 
arises because the relative distance between the charges depends on both their bouquet 
positions on the centerline and the local rotation of the bouquet about the centerline. 
To find this orientation one needs to integrate auxiliary equations of motion at each 
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step {Darboux^s vector), so the solution must be known before the relative distance 
between the charges can be computed. However, the solution explicitly depends on the 
relative distance between the charges, and so the closure of the system is problematic. 
These complications make an explicit derivation of equations of motion from Kirchhoff 's 
approach difficult, if not impossible. On the other hand, the approach we suggest here 
provides a very straightforward derivation. One need not deal with vectors, forces and 
torques in non-inertial frames of references, in writing conservation laws, and so forth. 
Instead, in the geometric approach, helices are treated in exactly the same fashion as 
straight lines, so stationary helical states and even their linear stability properties may 
be considered conveniently. 



2. Derivation of equations in SE{3) coordinates 

In this section, we derive the full equations of motion for a self-interacting rod in the 
discrete and continuous cases. This derivation is based on the exact geometric rod theory 
[23], that derives equations equivalent to Kirchhoff 's equations for elastic rods using 
symmetry-reduced variables. Since the derivation of the exact geometric rod equations 
for purely elastic rods is well-established, we shall concentrate on the corresponding 
derivation of the nonlocal equations in the language of the Lie group SE{3) - the 
Special Euclidean group of orthogonal rotations and translations in M.^. This familiar 
Lie group is a semidirect product of 5*0(3) and with the following definition of group 
multiplication: 

(Ai, ri) ■ (A2, ra) = (A1A2, Aira + r^) (1) 

for Ai,A2 G 5*0(3) and ri,r2 G M^. We use the notation 5e(3) for the corresponding 
Lie algebra. For the reader's convenience, we summarize the properties of this group, 



its adjoint and coadjoint actions Ad, Ad*, ad, ad* in [Appendix A[ More details can be 
found in, for example. 



Definition 2.1 A bouquet is a rigid, non-deformable assembly of spheres, each 
characterized by an interaction radius (for short-range interactions like Lennard-J ones) , 
and an electrostatic charge. 

A molecule, as considered in this paper, consists of a rod that represents the 
molecule's central axis, with bouquets attached at fixed base points along the rod. Each 
bouquet is characterized by the spatial position coordinate of its base r(s) G and 
its orientation A(s) G 5*0(3), which together define an element a{s) = (A(s),r(s)) G 
SE{3) = 50(3) X M^. We are interested in describing the effects of nonlocal interactions, 
which in this paper denotes all inter-atomic forces that affect parts of the molecule not 
immediately adjacent to each other. These are all the forces that do not come from 
elastic deformation, and include, for example, electrostatic and Lennard- Jones forces. 
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(4) 



To describe nonlocal interactions among these bouquets of charges, we define their 
relative orientation and position variables by the SE{'i) product 

= a~\s)a{s') = {A-\s)A{s') , A^\s){r{s') - r{s)) 

■= {^{s,s') , k{s,s')) , 

where ^(s, s') G S0{3) and k{s, s') E M'^. The position of the k-th charge in the bouquet 
whose base is at s along the center-line of the rod is given by 

Ckis) = ris) + Ais)r,,is) . (3) 

Here ^7^,(5) denotes a vector from the base point at s on the rod to the k-th off-axis 
charge of the bouquet as measured in a rigid Cartesian frame with orientation A(s). 
Then, the distance dkm.{s,s') between the k-th charge at s and m-th charge at s' is 
given by 

dkmis, S') = |Cfc(s) - Cm{s')\ 

= \r{s)+A{s)ij,is)-r{s')-A{s')rj^is')\ 

= \A-\s){r{s) - r{s')) + rj,{s) - A-\s)A{s')rj^{s')\ 

= \i^{s, s') + r]k{s) - ^(s, s')r]„Xs')\ = dkmC^is, s')) . 

Remark 2.2 In finding particular helical solutions and analyzing their stability, we 
consider the case of polymers where all the bouquets are exactly the same, so r]f^{s) and 
r]^{s') are independent of s and s' . Nonetheless, the equations of motion we derive 
in this section, (22) and (23), would also be valid for arbitrary dependence of the 
bouquet's geometry on position. However, the linear stability analysis performed later 
will explicitly use the fact that all bouquets are identical, and thus cannot be applied to 
the conformations of molecules with varying geometry. 

Consequently, the total energy of the rod due to nonlocal interactions among the 
charges that compose it is obtained in the continuous 

^ = Y.j U{dkm{^{s, s'))) ds ds' , (5) 

k,m 

for interaction potential U between individual pairs of charges. In the discrete case, the 
integral becomes a sum and s, s' are discrete indices, so that 

E= Yl U{dU^is,s'))). (6) 

s,s' ,k,m 



The corresponding contribution of nonlocal interactions to the action in Hamilton's 
principle for the dynamics of the rod is thus (in the discrete case) 



5. 



-Jj2 U{dU^{s,s')))dt. (7) 

s,s' ,k,m 

In order to compute the contribution of the mutual charge interactions to the dynamics 
of the rod one must take the variation of this nonlocal part of the action with respect 
to the charge conformation S. Upon denoting 

u{s) = a'\s)Sa{s) e st{3) , (8) 
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where 5c(3) is the Lie algebra of the Lie group 5* £"(3), we find 

s') = 5{a-\s)a{s')) = -u{s)E{s, s') + s')u{s') . (9) 

We now define a pairing (-, ■)tse{3) between the elements of TSE{3) and TSE{3)* - the 
tangent and cotangent spaces of SE{3) - as follows. If a vector U is tangent to SE{3) 
at point a E SE{3), and W co-tangent to SE{3) at the same point a, then a~^U brings 
the vector U to the identity element of the group, so a^^U G sc(3) is in the Lie algebra. 
Similarly, a~^W G 5c(3)*, so a scalar product between these elements can be taken as 
defined in the appendix. Thus, we define 

{U,W)^SE(3) = {^''U,a-'W)^^^,y (10) 

In what follows, we drop the subscript TSE{3) from the scalar product as not to burden 
the notation unnecessarily. 

In order to use the minimal action principle and derive the equations of motion, 
we need to take the variations with respect to S. The variation 5S defined in ^ is 
an element of TSE{3), as SE is tangent to SE{3) at the point S. Then, the derivative 
SU/SE is an element of TSE{3)* as it is co-tangent to the group at the same point E. 
The pairing 



can be now defined according to (10). Consequently, the variation of the nonlocal part 
of the action in ([T]) is given by 



6E 

s,s'k,m 



s,s'k,rn 



:ir 



'■,s' ,k,m 

in which the last step uses the relation E{s', s) = E~^{s, s') obtained from the definition 
of S in (g. We shall now proceed with the computation of the variations of action with 
respect to all dynamical quantities and thereby obtain the equations of motion. 



Velocities. The local part of the Lagrangian in Hamilton's principle for the dynamics 
of the rod is written by introducing the left-invariant variables 

= a-^a = (A , rf^ (a , = (^A-^A , A-^) G se(3) , (12) 

as velocities taking values in the Lie algebra se(3). 



Helical states of nonlocally interacting molecules 7 

Elastic deformations. In the continuous case, the invariant variables that describe 
elastic deformations are, 

A = a-V'Gse(3). (13) 

In the discrete case, following the Moser-Veselov method for numerical discretization of 
rigid body dynamics [2S] we set, 

A = (7-i(s)a(s + l) e 5^(3), (14) 

where s = 1,2,... is the discrete index labeling a given base. The elastic part of the 
Lagrangian will then depend on A. 



Compatibility conditions. The compatibility conditions are obtained, in the continuous 
case, by using equality of cross derivatives, so that cTst = Cts- Differentiating (12) with 
respect to s and ( 13 ) with respect to t, then subtracting yields 



-Xfi + fiX := [fi , X]s,(3) = ad^X , (15) 



9/i dX 
ds dt 

where [fi , X]sc(3) is the commutator of fi and A in 5e(3). In the discrete case, we write 
(14) as a{s + l,t) = cr(s, t)A(s, t). Differentiating this condition with respect to time, 
we get 

dX 

— (s) = A(s)/i(s + l)-M^)A(s), (16) 
or in terms of se(3)-algebra quantities only 

x-'—{s)=^l{s + l)-Ad,-^^l{s). (17) 

Dynamical equations: variations In the continuous case, the variations of the velocities 
on the Lie algebra se(3) satisfy, 

du 

6fi = -ufi + fiu + — . (18) 
In the discrete case, this is replaced by a variation on the Lie group SE{3), 

X-^6X{s) = -Adx-iiy{s) + iy{s + l). (19) 
For the discrete case, one may then compute the variation of the Lagrangian as. 



sx 



5Z (^"'(^) aY^^) ' -^^x-HsAs) + u{s + 1)) (20) 
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In continuous time, the Euler-Poincare equations emerge from the following direct 
computation as in [211 122] . 





5/.'"^-'^+9^r^' = y J \[-dt^<)r^'''r'^'- 

Likewise, the equations of motion in continuous time for the spatially discrete nonlocal 
rod, may be written in SE{3) coordinates as 

| + ad;)^-Ad.-,„(A-w|(,))+A(.-l)|^(.-l) 

s',m,fc 

In the continuous similar calculation gives (note that now A = cr ^a' G se(3)): 

-E /-^(■'■s')H-'(s,.5')+H(s,.5')^(s',s)d/ = 0. 



Remark Notice that (23) are exactly the equations derived earlier in [211 122] • However, 
the discrete equations (22) for nonlocally interacting molecules are new, as far as we 
know. 



Lemma 2.3 For an arbitrary 'potential U{d), equations (22) and (23) reduce to 
algebraic equations for helical solutions. The solutions of these algebraic equations are 
stationary helical shapes. 

Proof Helical configurations in the discrete case are obtained from taking a given 
element a G SE{3) and defining a{s) = a^. In the continuous case, a given element 
A G se(3) generates a{s) from the differential equation cr'(s) = o"(s)Ao (exponential 
map). In either case, 

S(s, s') = a-\s)a{s') = a{s' - s) . (24) 

Then, 

dkm{s, s') = dkmC^is, s')) = dkm{cr{s' - s)) , 

and therefore 

where Di is some TS'-E'(3)*-valued function that depends only on the difference between 
s and s'. Thus, the right-hand sides under the sum (22) and in the integral in (23) 
depend only on the difference between s and s', and are therefore constant. Likewise, 
the left-hand sides of these equations are also constant, as = and A = a G 5-^(3) in 
the discrete case and A = cxo G 5c(3) in the continuous case. Thus, integro-differential 
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equations (22) and (23) reduce to algebraic equations for helical solutions, as long as the 
potential between two charges depends (in an arbitrary fashion) only on the Euclidean 
distance between those charges. □ 

Of course, this result had been well-known by molecular biologists for at least 
50 years; in their famous paper Pauling, Corey and Branson noted that ". . . is 
likely that these [helical] configurations constitute an important pari of the structure of 
both fibrous and globular proteins, as well as synthetic polypeptides'". As it turned out, 
the helical structures are amazingly robust; they appear for large variety of molecules 



with widely ranging elastic properties and charge distributions. Lemma 2^ provides a 
mathematical model for this well-established empirical fact. It is worth noting however 
that the molecules that are locally helical almost never stay in the shape of a perfect 
straight helix if they are sufficiently long. This folding of molecules is driven by a 
combination of the elastic and non-local forces and plays a very important role in 
the ensuing functionality of molecules. The instabilities of the helical states will be 
considered later using geometric approach. 



Remark 2.4 Lemma 2^ includes, as a particular case, results from previous work on 
the helical solutions of the Kirchhoff's rod (no nonlocal terms) or uniformly charged 
Kirchhoff rod JT^ . with the charges'^ positioned at the rod's axis. In that case, the rigid 
charge conformations are positioned at r/(s) = 0, and thus 

d{s, s') = — s)| . 

In this case, the distance d is independent of the mutual orientation ^{s,s'), which is 
the first part of SE(3) group element a~^{s) a {s') := (^(s, s'), k;(s, s')). That is, charges 
arranged on the axis produce no torque. 

With this particular simplification, the computation still proceeds in the same way 
as before and reduction to algebraic equations still holds. Of course, more work is needed 
for each particular potential to demonstrate that the algebraic equations actually have a 
solution, especially if complicated elasticity laws in the rod are assumed. Thus, Lemma 



2.3 provides a simple justification for helical solutions that have been ubiquitous in the 
previous literature. However, our result goes further: helical shapes also allow one to 
search for helical solutions as solutions of algebraic equations for arbitrary charges off 
axis and for a charge distribution at each point of the molecule, possibly generating 
torque in the molecule. 



3. Application to a linear rod with naturally straight conformation 

This section demonstrates how to find helical stationary shapes for charge bouquets 
positioned along an infinitely long, naturally straight molecule. 

To illustrate the general ideas described above, we shall consider a molecule that 
consists of a chain of centerline atoms connected with elastic springs, and charges 

^ Here, we use the word "charge" somewhat loosely to denote arbitrary potential interaction, like 
Morse or Lennard- Jones, and not just electrostatic potential. 
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that are attached to each unit of the elastic chain with a rigid charge bouquet. In 
the undisturbed (base) configuration, the molecule is straight. In what follows, we 
select a charge bouquet containing two charges of qi = ±0.3e, where e is the charge 
of electron. This arrangement approximates any molecule which has a constant dipole 
moment perpendicular to the axis, for example, vinylidene fiuoride oligomers (VDF) 
[27] . More general charge configurations such as quadrupoles etc. can be incorporated 
by considering more general bouquets. 




Helical Symmetry 



Figure 1. Top left: an example of a charge bouquet with four atoms. The electrically 
charged atoms are shown in red and green (positive or negative). Grey atoms are 



neutral, but still interact with other atoms through Lennard- Jones potential (34). The 
molecules we consider are long chains of these bouquets. 



We assume electrical charges exhibit a screened electrostatic interaction 



Eij,c{s,s' 



Qi[s)qj[s 



-dij{s,s')/\ 



i-j\S, s 



(25) 



where dij is the distance between charges qi at s and qj at s', and A denotes the Debye 
screening length. In what follows, we consider values of A corresponding to various ionic 
strengths in the solution, from / = 0.001 M/1 (close to de-ionized water) to / = 10 M/1 
(an order of magnitude higher than sea water). Debye length A and ionic strength / are 
related by 



A 



eoerksT 



(26) 
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where eo is the permittivity of space, is dielectric constant of the media, is 
Boltzmann's constant and Na is Avogadro's number. 

We seek stationary states that are invariant with respect to an affine (hehcal) 
transformation of space, which rotates a coordinate frame F by a matrix A, and performs 
a translation of an arbitrary vector r as 

(F, r) ^ (AF, Ar + a) . (27) 

In particular, we seek solutions with bouquets spaced uniformly on a helix of radius R 
and pitch C, with the rotation being parallel to the (x, y) plane. By pitch, we mean 
that after making one period of rotation, the atoms are moved by the distance C along 
the 2;-axis. The helix is then defined parametrically by 

x = Rcost, y = Rsint , z = Ct , (28) 

or explicitly 

x = Rcosz/C, y = Rsmz/C, (29) 

where R> and C ^ 0. The helix is right-handed for C > and left-handed for C < 0. 
We use R and C as parameters in our calculations and visualizations. A continuous 
helix is a curve that is invariant under a set of transformations consisting of a rotation 
A about the z axis and a translation a, parameterized by z, 



K{z) 



cos(^) sin(§) 
-sin(§) cos(§) 
1 



a U 



i?sin(§) 

Rcos (^) 
z 



(30) 



For a discrete helix, we consider values of t, (respectively, z) that are discrete multiples 
of some value T (respectively, z/C)\ 



t = nT 



resp. 



n 



C 



where n G Z . 



(31) 



Bouquets on the helix are oriented such that the helical invariance described above 
is maintained. Once a single bouquet is specified, at say 2; = 0, (30) generates the entire 
helix structure. Let A G 50(3) denote the orientation of the initial bouquet at 2 = 0. 
In what follows, we limit A to a twist about the tangent to the helix at z = 0, as shown 
in Fig. [2| allowing a single angle a to characterize bouquet orientation. 

-C ■ ^, 2-kR 



A{R,C,a) 



cos a 

c 



sma 



-2ttR 



sma 



sma 



iw^R'^+C^ cos a 



47r2_R2+C2 
2ttRC r 



V47r2i?2+C2 " 

27rRC q_cosa) 



47r2i?2+c'2 ' 



COS a 



47r2_R2+C2 

47r2R2 cosa+C2 
47r2ij2+c'2 



(32) 



The entire helical structure may then be generated by repeated application of (30). 



4. Energy of a helical conformation 

This section explains the interaction energies of a realistic molecule to be used in 
calculations, both the local (elastic) parts and non-local (electrostatic and Lennard- 
Jones). As discussed above, we investigate the particular example of a molecule of 
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Figure 2. Effect of twist transformation - the twisted bouquet is highlighted. 



bouquets with two opposite charges ±g positioned on either side of the polymer axis, 
each a distance I from that axis. This simphfied model describes a polymer with a 
constant polarization perpendicular to the polymer's axis. The interaction energies are 
defined as follows. 

Elastic energy. We consider elastic energy as quadratic in bond bend angle, 



where /i is a spring constant, and A0 is the bond bend angle with respect to its 
unstressed orientation in the intrinsic frame, assumed straight. For naturally helical 
molecules, A0 would be the value of rotation angle for a bouquet about its rotation 
axis in the intrinsic frame from its equilibrium value. In terms of global geometry 
that would mean that bend energy is punishing creation of local curvature. A value 
of /X = 3.025 ■ lQ^J/{ra(P)Mol = 5.022 ■ IQ-^^ J/ra(P is typical for bend rigidity for a 
carbon-carbon bond. The bend depends (through a complex algebraic formula that we 
do not present here) on the radius and pitch (axial distance traversed on each rotation) 
of the helix. 

Lennard- Jones energy. To prevent self-intersection, we introduce a truncated Lennard- 
Jones interaction between charge bouquet nodes with equilibrium distance d^ and 
potential well depth e. When the centers of two bouquet nodes are a distance d from 
each other, they experience a potential given by 



'elastic 



(33) 




(34) 
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Charge potential energy. As mentioned when the model was introduced, each bouquet 
has charges of qi^2{s) = ±0.3e that interact with each other through the screened 



Coulomb potential (25) with Debye length A. 



Total energy of a given conformation. Given a bouquet configuration and values for R, 
C, and a, we compute the total energy by choosing a reference bouquet, then working 
outward along the helix in both directions by helical symmetry. For each bouquet we 
reach, its contribution is the sum of the Lennard- Jones and Coulomb interactions of 
the bouquet nodes with the reference bouquet nodes. The sum of these, plus the elastic 
energy corresponding to the given R and C give the energy per bouquet (energy density) 
of the given conformation. 

As an example, we demonstrate the particular energy landscape for / = O.OOIM// 
in Fig. |3j The vertical axis is helix radius R, the horizontal axis is helix pitch C, 
with {R, C) = (0, 0) in the upper left corner. The energy scale bar is shown on the 
right, with lowest energies in black, highest in white. The left edge and the upper left- 
hand corner contain most of the energy minima. The minima along the left-hand edge, 
corresponding to small pitch C (on the order of the size of the bouquet), corresponds to 
the helices solutions that are being "ratcheted" into tighter and tighter conformations 
as we proceed upward (smaller R). The values of the energy in this graph are minima 
over all twists (values of a) of the bouquet, as shown in Fig. [Tj The local minima of 
this energy landscape are stable helical conformations. 

In order to further elucidate the structure of the energy landscape. Fig. |4] presents 
a blow-up of the two boxes in Fig. [3] Yellow dots correspond to energy minima and 



thus show stationary conformations. For completeness, in Appendix B we show the 
change of energy landscape as the ionic strength of the surrounding solution varies from 
0.001 M/1 to 0.016 M/1. Since we assume that the base conformation of our molecule 
is linear, a weakly ionized solution leads to longer electrostatic interaction, resulting 
in richer structure. Helical conformations become increasingly unlikely for strongly 
ionized solutions that make long-range electrostatic forces weaker. For strongly ionized 
solutions, elastic forces tend to unwind the helices into straight lines. 

As shown in Fig. |4| even this simple molecule, when deformed into helical 
configurations, shows quite a complex and intriguing energy landscape. While the 
details of the landscape depend on the individual molecular parameters, the presence of 
"ratcheting" helical states is, we believe, typical for all molecules of this type. In this 
figure, several examples of the helical conformations corresponding to different energy 
minima are presented as inserts. There are two types of helical conformations. The 
"twist" type is obtained by twisting a base straight conformation, and is characterized 
by small radius R and increasing pitch C. In a given energy landscape, there are only 
finitely many conformations of this type, as increasing the twist beyond certain level 
will cause elastic forces that are too large to be balanced by electrostatic attraction. We 
show two examples of such twist conformations in the top part of this figure. Another 
type of helical solution is given by the ratcheting states, characterized by a pitch on the 
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Figure 3. The energy landscape for a given ionic strength / = O.OOlAf/L The 
vertical coordinate is radius R and the horizontal coordinate is pitch C . Energy is 
also dependent on the rotation of the bouquet around the axis a, but this coordinate 
has been projected out onto two dimensions by taking a minimum over all rotation 
angles for the bouquet. That projection preserves the minima of the energy. Notice 
that this energy plot provides a complete analysis for all helical conformations. Helical 
conformations that are potential minima are concentrated on the left-hand side (small 
pitch, increasing radius) and upper left-hand corner (small radius, increasing pitch. 
Below, Figure [4] shows blow-ups of the framed rectangular regions. 

order of bouquet's size and increasing radius. This is shown both in the top and bottom 
parts of the figure. There are infinitely many of these ratcheting states, as the radius can 
increase (in principle) to arbitrarily large values. Energy minima then occur every time 
the opposite charges line up on sequential rolls of the helix. That explains the regularity 
of occurrence of those minima. They occur every time the helix circumference increases 
by the distance separating the bouquets, /. This leads to an increase of radius between 
neighboring ratcheting conformations as Ai? = //(27r). 

To illustrate how varying ionic strength of the solution changes the energy landscape 
of helical conformations, we compute 10 energy landscapes varying from / = O.OOIM/Z 
to / = 0.16M/L These energy landscapes are presented, in an animated fashion and 
with the same orientation and axes as in Fig. |3} in the supplement. From this, it is clear 
that the variation of ionic strength and hence Debye length, has a profound effect on the 
energy landscape. As we see below, the ionic strength is also of particular importance 
for linear stability. 

5. Multiple helical conformations 

In this section, we explore generalizations of the helical shapes considered above. We 
show how to compute a particular molecular shape that we refer to as a 2-helix (and, 
in general, n-helix). By 2-helix we mean a conformation where a group consisting of 
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Figure 4. Blow-up of energy landscape for the boxes in Fig. [3] The energy scale has 
been adjusted to show more detail, and is shown on the right with lowest energies in 
black, highest energies in white. The vertical coordinate is radius R and horizontal 
coordinate is pitch C . Top: upper left corner box in Fig.|3] Bottom: box in the middle 
of left edge in Fig. [3j Several examples of the helical conformations corresponding 
to different energy minima are presented as inserts with yellow arrows indicating the 
corresponding energy minima. 

two bouquets with some relative position and orientation is repeated along a helix. 
Alternatively, a 2-helix can be seen as two helices of single bouquets arranged so 
the polymer axis passes through bouquets from alternating helices over its length. 
The relative orientation of bouquets in the group can be arbitrary. Such structures 
have been suggested and analyzed in |28] in the context of describing given molecular 
configurations: /3-helix structure (2-helix) or collagen (3-helix). These structures may 
be obtained as natural energy minima for a given molecule. In particular, we present 
an example of an equilibrium configuration of a 2-helix using the bouquet considered in 
the previous section. An extension of these ideas is possible for n-helices with n > 2. In 
general, ra-helices are combinations of n helical conformations, with 2, 3, . . . , n being 
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obtained from the first one by a shift, rotation and a twist of the base bouquet. A 
detailed map of the energy landscape for an n-helix is problematic even for moderate n, 
because the number of parameters involved in defining the conformation of such a shape 
may increase prohibitively with n, as the position and orientation of each bouquet is 
required in defining a repeated group in an n-helix. 

For a 2-helix, we optimize the energy over both the shape of the helix and the 
relative configuration of two bouquets. A minimal energy over these variables will, by 
the helical symmetry, provide an exact stationary conformation for the whole molecule. 
Two examples of such conformations are given in Fig. |5j Details of calculations of more 
general multi-helices will be discussed elsewhere. 




Figure 5. 2-helix conformations for the ionic strength 1= 0.001 M/1, and red/green 
charges in each bouquet of ±0.8e. The elastic centerline (not shown) goes through 
gray spheres that are the centers of the bouquets. The conformation is not a perfect 
helix, but instead consists of two helices, in which the second helix has undergone a 
twist, shift and rotation with respect to the first one. The distance between the bases 
of the bouquets (gray spheres) is 2 A. 



6. Linear stability analysis 

This section derives the dispersion relations for the linear stability of helical polymers. 



based on the linearizations of (22) and (23) about a helical state that is assumed to be 



a stationary solution. As we show in this section, the use of geometric methods yields 
an exact dispersion relation. 

The energy landscape for charged polymers with nonlocal interactions is complex, 
even when limited to helical shapes, because of the large number of energy minima 
in that "helical universe". It is natural to pose the question of linear stability and 
selection of the stationary shapes. Here, the geometric approach allows us to find the 
exact dispersion relation for the stability of the helical shapes. As far as we know, no 
such work has been undertaken before, perhaps because of the complexity of the linear 
stability analysis when using traditional methods. 

By using an SE{3) method that naturally accounts for the helical geometry, one 
may investigate the linear stability of a general helical polymer with torque-inducing 
non-local self-interactions. Studies of the stability of helical elastic rods have been 
undertaken before, with most studies concentrating on the linearization of the traditional 
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Kirchhoff equations about the stationary helical states [121 EHl EQl [31], |32l [33]. The 
focus of these works was on the stability analysis based on increasingly complex elastic 
properties of the rod. Alternatively, the work [M] investigated the stability of elastic 
rods using the exact geometric rod theory and applied it to DNA dynamics. All these 
works have used the continuous model of elastic rods. Our results differ in two ways 
from previous studies. First, our results are formulated for spatially discrete rods, 
as discussed in the previous sections of this paper. Second, and more importantly, 
our stability analysis includes non-local interaction of charges that in general occupy 
positions off the axis of the elastic rod. The forces on these charges generate torques 
acting on the rod's centerline. These torques are absent when one considers purely 
elastic rods, or when the charges are distributed only along the centerline. 

This section shows that the presence of the torques due to non-local interaction of 
off-axis charges generates an instability of the rod and that the instability appears even 
for the simplest possible states - the linear rod. The instability due to nonlocal torques 
is new, as far as we are aware. 

Obtaining similar results using Kirchhoff's rod equations would be problematic. 
The main difficulty consists in finding the Euclidean distance between two arbitrary 
points on a rod's axis using coordinates intrinsic to the rod. In addition, the present 
computation shows the geometric origin of the exact dispersion relation for arbitrary 
helical configuration of a rod (without charges), something that was noticed already in 
|12j . In our opinion, the present method based on exact geometric theory is more 
straightforward, algorithmic and compact than the corresponding linear analysis of 
Kirchhoff equations. 

Suppose we have a helical configuration arising from successive repetition of a given 
element a G S'i?(3), so that a{s) = and S(s, s') = a'^''^. Let us linearize about this 
solution, so 



As it turns out, this substitution leads to an exact dispersion relation for ai. We define 
the 0(e) perturbations in velocity and deformation as follows. First, perturbations of 
the velocity are determined, as 




(35) 



The linearization written symbolically as a = ctq + is defined by 




(36) 





(37) 



In the discrete case. 




(38) 



where 



Ai = -o-g ^(s)cri(s)a + cro ^(s)(Ti(s + 1) 



^i(s)a + z^i(s + l) Gse(3).(39) 
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S cr S 



r and finds, 



(40) 



In the corresponding continuous case, one denotes Aq 

A(s,t) = r + eAi + 0(e2) 

Ai = -ffo^Vir + crQ^a[{s) = -ij^T + i;[ e se(3) . 

For simplicity (in order to keep the formulas compact), we shall assume for the discrete 
case that 

1^ = Ho + e//ii + . . . , X~^^ = Ko + eK,{s) + .... (41) 

O/i OA 

Next, we compute the linearization of the nonlinear terms due to elasticity. In order to 
find the linearization of the sum of terms 



Adl-i(,)(A-i 



(42) 



in equation (22) we utilize the following proposition. 



Proposition 6.1 Suppose G is a Lie group with Lie algebra q and (-,■): g* x g M 

is a pairing between the Lie algebra and its dual. Suppose A{e) E G is a curve in G with 
A{0) = Aq and Aq^A'{0) = a E q , a{e) G Q* , and a'{0) = ^. Then (see, for example, 

a p-60) 

d 



Ad^^i (^^ - ad\ao) . 



(43) 



Note that X~^6l/6X has the physical meaning of the local stress in the body coordinate 
frame at point s. 
Then, writing 



where , Ki G sc(3)*. Hence, the linearization of (42) is computed as follows 



d_ 
de 



e=0 



-M\-r [Ko - ad;^iri(s)) {s) + K^{s-l), 



(44) 



where '4^i{s) := Ag ^(s)Ai(s) =G 5e(3). 

The linearization of the nonlocal terms in equation (22) is less straightforward and 
will be outlined in its own section below. For now, we assume it is possible to compute 
that linearization, and it is described by some linear operator L(S(s, s'))ipi{s), which is 
defined as follows. Consider an arbitrary rj G se(3) and define the scalar function of s 
by the following pairing. 



Us) 



^(s,s')S(s,s') + S ^(s,s')^(s',s) 



vi)ds' . 



(45) 



The nonlocal term takes values in the space se(3)*, so the pairing in (45) indeed defines 
a scalar function. We need to find its linearization with respect to ip = a~^5(y G se(3). 
For this, we compute the derivative of I{s) with respect to S according to 



51 



(46) 
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calculation, we need to express 



where ^/'i = S ^Si G se(3) is the linearization with respect to S. To complete this 



(47) 



Li in terms of ipi. This step proceeds as follows. 
The linearization of s') in (24) gives 

= a'^'"^ + eHi(s,s') + 0(e2), 

where 

Si(s,s') = - a"''cri(s)a^'"'' + a~Vi(s') 
= -V^i(s)a^'-^ + a^'->i(s') 
= a^'-^(V'i(s')-Ad,_.Vi(s)). 
Consequently, the quantity S^-'^Si is given by 



and we have 
51=1 - Ad* 



-Ad^s-.'V^i(s) + V'i(s') 



5S 



, , 5I{E{s',s),r]) 
+ ^{s , s) — , Ipi 



(48) 
(49) 

(50) 



Finally, since I{s) in (45) is a linear function of an arbitrary rj, re-arranging expression 
(50) into a scalar product of an se(3)*-valued function with t] will give the desired 
linearization operator L(H(s, s'))'?/'i(s), from 

51 =: {UE{s,s'))Ms),v)- (51) 
The equations simplify further upon noticing that for a stationary helical solution, 
a = aoa'^ where a G SE{3) is a given element, so the expression Aq = aQ^{s)a{s + 1) = a 



is independent of s. Then, the linearization of equation (22) in the discrete case is 
,M - Ad:-x (iri(s) - ad;^is)Ko) (s) + K,is- 1) 
= 5^L(H(.,.'))'0i(3). 



(52) 



s'm,k 



It is natural to posit the following ansatz: 

Ki{s) = [j(^i(s + 1) - V^i(s))] a' := M^s {j{i^i{s + 1) - Vi(s))) , (53) 

where ipi{s) G se(3) and J : 5c(3) -> 5e(3)* is a linear operator having the physical 
meaning of the rigidity matrix. Notice that the linearized system of coordinates is 
written at the point s, but it encounters the value of the stress at the point s — 1. In 
order to connect this stress with the coordinate system at the point s, we will need to 
transform the coordinates to s — 1, by shifting one step forward on the helix. We thus 
need to compute Ad*-i of the term evaluated at s — 1, i.e. 



J(^i(s)-V^i(s-1)) 



Then, the linearization of the discrete case gives 
92 



li^i - Ml 



j(^i(s + l)-2^i(s)+V^i(s-l)) 



(54) 



s'm,k 
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In the continuous case, the corresponding hnearization of equation (23) gives 

-^M+(-- + ad*)j^; + ad;^r = 5^ / L(S(s,s'))V^i(s)ds'. (55) 

m,fc 

Here, again i\} = Aq "'^(s)Ai(s) G 5c(3). 

Further simphfication can be obtained for the nonlocal term for the stationary 
helical state a{s) = aoa'^. The invariant variable S = a~^{s)a(s') = a^~^' depends 
only on the difference between s and s'. Thus, all of the derivatives of the potential 
energy with respect to S when evaluated at the helical configuration depend only on 
the difference between s and s'. In other words, we have 

with e se(3)* for any a e 5c(3), so D2 : sc(3) sc(3)*. 

Since all the functions on the right-hand side depend only on the difference s — s', 
the integrals or sums become convolution integrals. Fourier transforming then allows 
the exact dispersion relation to be obtained, as follows. Let us consider 

^1 {s,t) = Se-''^'^'^' , S G se(3) . (57) 

Here s is an integer and k is the dimensionless wave number, measured in the units of 
2tt/Iq, where Iq is the distance between the elements of the helical chain. Consequently, 
assuming S is real, the linearized equation (54) gives the following dispersion relation: 

u'lS - 4( sin^ ^)Ad:-. ( - ^d*sKo) = Us')S . (58) 

s' ,m,k 

Note that in the absence of non-local interactions, the basic helix must be unstressed, 
Kq = 0, so = a are given by the generalized eigenvalues of the problem 

4 sin^ ^Ad;;.! (^JS^ = IXS . (59) 

From physical principles, we require all the generalized eigenvalues of the matrices J and 
/ to satisfy A = > 0, so that all purely elastic helices in stationary conformations 
are neutrally stable. The spatially discrete dispersion relation (58) converges to the 
dispersion relation for the continuum case in the limit A; — 0, a — )■ Id5£;(3), which is 

fc^Ad^-i (js) = IXS . (60) 



Again, the right-hand side of ( 58 ) is a function of s' — s only, while the left-hand 
side is a constant. Upon summation over s', the dependence on s disappears and the 
dispersion relation is obtained by setting s = on the right-hand side: 

^^IS - ( sin^ 'l)AdU {jS - ad^i^o) = Us')S . (51) 

s' ,m,k 

The right-hand side is a linear operator acting on S. Instability corresponds to 
generalized eigenvalues u of equation (61) having a positive imaginary part. As it 
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turns out, all eigenvalues A = are real, so it is enough to identify the case A < as 
instability. However, the consideration of discrete rods puts an interesting spin on this 
problem that we will consider below. 



7. Computation of derivatives for potential energy 



Equation (61) provides the stability analysis for an arbitrary interaction potential 



U{d). However, for particular applications it is important to explicitly compute 



the linearization operators in (61). Again, geometric methods will be advantageous. 
Therefore, this section computes the linearization of the nonlocal terms for the dispersion 
relation. We believe it is advantageous to show this computation in some detail, as it is 
not trivial. 

We shall perform the computation only for the discrete case. The continuous case 
is derived similarly with the change of sums with respect to s and s' into integrals where 
necessary. One has to be careful here, as we need to take derivatives of quantities that 
take values in TeSE{3) and TeSE{3)*. The most straightforward way, least likely to lead 
to a mistake, is to define a corresponding scalar functional by bringing these quantities 
to the Lie algebra and then pairing them with the corresponding fixed element from 
the dual. The derivatives will then be given by whatever term is paired the chosen 
fixed element. This is akin to the weak computations of functional derivatives, only 
performed with geometric quantities. First, notice that 

^(., s')E-\s, s') = Ad|-.(,,,,) {e'\s, .')^(., .')) , (62) 

and 



.s, s 



(63) 



Thus, we start the computation of the linearization operator D2 for the nonlocal 
term by linearizing the expression G se(3)*. Let us consider S = (^, k) with 

= {C,~^, an arbitrary element (0, ^) G se(3) and a scalar functional 

h = [r. — — — ,{ti,a)J. (64) 
The linearization operator D2 is therefore defined as 

2^2(0,*), (m,«)) := (S-I^^i^, (0,*)). (65) 

Here, we introduced the natural pairing between two elements (a, b) G TeSE{3)* and 
(«,/3) GT,^E(3): 

{{a,b), {a,(3)) = ^tT{a^a) + b- (3. (66) 

Given 

4m(2(s,s')) = dkm{^,K) = \k + r]i^{s) ~ ^{s,s')r]^{s') \ , 
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U'idkm) 



tr(| - ^dkm <H) rij'^fij + ^dkm) ■ a 



(67) 



km 



since {fi)ij — e^jki^k by the definition of the hat map. Here, we have defined 

U'{dkm) 



d 



km 



(68) 



/2(S, /U, a) := -{i ^dkm X ImV • M + ^^km) ' OL . 

In order to compute the hnearization D2, we proceed as follows. For (</>,*) G 5e(3), 



calculate 



Q'{d 



km, I 



d, 



km 



/2(S,/i,a)/2(S,(/),*) 



(69) 



We still need to compute the variational derivative of hi'^, l^-, ex). The only part of I2 
depending on S = (^, k.) is the quantity 



Then, 



r'^r'('*+%)-a,z/' 

X ^ 1 

= --tr < r'^('« + %) • e« , ^) - 2t^(« <^ + ' ^) 



(a X (/t + 77fe)) • * = ((k + r7fc) X *) • a . 



Similarly, 



r'^r'(/* + r?,) • a , ^) = (r'^('* + r/,) • , ^) = « • * ■ (71) 

The derivatives of ^~^(k. + r/^) x t)^ are computed similarly using standard properties 
of vector cross products. For brevity, we shall not present these calculations here. The 
final answer for 1)2,^ is given by collecting the terms proportional to /it, and coefficient 
of D2,a is given by the terms proportional to a. The operator D2 is thus given by 



^2,«(*,</>) 



d 



km 



-/2(5, </>,*) (77m xr'dfc„) 



+ Qidkm)Vm X ( - r'(« + %) X (/) + , 

Q'{dkm) J 



(72) 



''km 



-/2(H,0,*)r'dfem 



+ Q{dkm) ( - C ^(k + r7fc) X + . 
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From (72) one notices a very interesting relationship, namely, 

Z}2,M(*»=^n.X^2,«(*». (73) 

The linearization operator is computed from D2 as follows (using S(s, s') = a^'~^ 
and iii{s) = Se^'"'): 

L(S)(0, *) = -Ad|_i (d2((0, *) - ad^^ „)Di(S)) (s, s') + (^2(0, *)) (s', s) , (74) 
where 

/^i(2) = H-^^ = _ (e-M,^ X rjj^ (rM.„)) . (75) 

It is also useful to outline the formula for the change of variables s ■h^ s' that forms the 
last term of the linearization operator L. Under this change, the operators Ad and Ad* 
change their form, and it is essential to perform this transformation correctly. Since 

(0,*)(s,s') = -AdH-i(.,.o^i(s) +^i(s') = -Ad^.-.^e^'^^^ + e^'^'^'^, 
an exchange of variables s ^ s' gives 

(0, *)(s', s) = -AdH-i(y,.)e''=^>i(s') + ^i(s) = -Ad„..-.e''=^'^ + e^^'S . 
Thus, the final expression for the linearization of the nonlocal term is 

L(s)V'i = -Ad* \D2ia''){ - Ad,_..5 + e*^^'^) 



^'^(-Ad _^,S+ei''s's)^'i-(^) 



8. Numerical stability of a linear polymer 



D2{a"''){-e"'''Ad^s'S + S). 



(76) 



In order to apply the general method of geometric linear stability of helical polymers 
derived in the previous section and show how our theory applies to a real-world case, we 
will solve the problem of linear stability in an example of a naturally straight, untwisted 
polymer in its unstressed configuration. The reason why we choose this polymer is the 
relative simplicity of the formulas, where all the Ad and Ad* operators are identities. 
Also, this is exactly the PVDF polymer considered in the Section |4| with the linear state 
being the most basic energy state of the molecule. It is thus interesting to determine 
the conditions for the linear state to become unstable, so other helical states computed 
in that section can be achieved. We shall note that the stability of the helical states 
computed in Section |4] can be considered analogously using the results of the previous 
section. Here, however, we shall avoid doing this as it will make the paper unnecessarily 
complex, due to the large number of helical states we have described. The stability 
of more complex stationary conformations, including the 2-helix, will be considered 
in forthcoming work. Thus, for the sake of simplicity, we shall only concentrate on 
the stability of a polymer that is perfectly straight in an unstressed configuration, as 
illustrated in Fig. [6j In this section, we also explain the difference between linear 
stability of a continuous and discrete polymer, and explain why a short enough polymer 
may exhibit stability whereas a long polymer will be unstable. 
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We assume the 2-charge bouquet with a charge of ±g on each end, resulting in 
a constant dipole moment perpendicular to the axis. The charges interact through 



a screened electrostatic interaction (25) and Lennard- Jones interactions (34). The 
charges are positioned away from the axis at the distance Ic = lA, and the distance 
between the charges is Iq = lA. We take electrostatic charges to be 0.17e, leading to the 
dipole moment of 5.33 x 10~^°C ■ m ^ 1.63D (Debye units), which is slightly smaller 
than a PVDF polymer having similar polarization structure. Note that this charge 
is smaller than the value q = 0.3e taken in the computation of the stationary states, 
since for q = 0.3 the subcritical bifurcation occurs for the values of ionic strength of 
about / ~ 200M/1 which is orders of magnitude larger than any values of I achievable 
experimentally. Thus, the linear polymer with charges q = ±0.3e at the ends of bouquets 
will be inherently unstable for all viable experimental conditions. 

8.1. Setup of the problem 



Figure 6. A particular example of charged rod with repeated configuration of two 
charges, plus (top) and minus (bottom) of g = 0.17e that interact through a screened 
electrostatic and Lennard- Jones potential. 



Remark 8.1 A it is well known that the stability of finite Kirchhoff rods strongly 
depends on the end conditions imposed on the edges of the rod [35], p. 81. The issue 
of choosing the right boundary conditions is a delicate one and, as far as we know, 
not entirely understood even for the Kirchhoff rods. Since the focus of this article is 
the investigation of the effects of nonlocal terms, we shall assume that the boundary 
conditions at the edges of the rod are such that Fourier transform analysis can be applied. 

For the purpose of this paper we shall assume the simplest possible shape of the 
6x6 elastic tensor J: we take J to be diagonal, with the values of 3 first diagonals 
= 3.025- 10 V/(rad2)Mo/ = 5.022- 10"^ V/rarf^ being the twist rigidity of a C-C bond, 
as outlined above, and the values of the last 3 diagonals (stretch rigidity along different 
directions) being /i/^. Such a choice of the elastic constants achieves J = Idexe in the 
dimensionless units. All lengths are then expressed in units of Iq. 

For the purpose of this paper, we take the inertia tensor to be 

Iq = diag{moll,moll,moll,mo,mo,mo), 
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where mo is the mass of the charged atom at the end of the rigid bouquet. In reahty, 
/ will be a symmetric positive-definite tensor depending on the exact nature of the 
polymer selected. Selecting the time scale r = a/TV"^ sets all the coefficients of the 
temporal and elastic terms exactly equal to unity. It is convenient to choose the unit of 
electrostatic charge as 



= a/ ineofilo ~ 1.514e . 

The value is chosen in such a way that two charges separated by /q interact with 
potential ±/i. The dimensionless Lennard- Jones amplitude is e//i ~ 2.90 x 10"'^. 

Our theory is also applicable for more complex values of elasticity tensors. However, 
the more complex elastic properties of the rod may themselves lead to instabilities, as 
earlier works show |2H1EI]- Thus, we shall assume the simplest possible elastic tensor in 
order to concentrate on the appearance of instabilities due to the long-range interactions. 

Limiting the considerations to unstressed linear polymers provides rather 
substantial simplifications in the expressions for the nonlocal terms. More precisely, 
the following simplifications hold: 

a = Parallel shift along the rod's axis by 1, 

Ada" = Identity in sc(3). 

Ad*. = Identity in se(3)*, (77) 

Kq =0 (no stress in the basic state), 

Di =0 (no twist in the basic state). 

Using this information, the dispersion relation u{k) can be now directly computed from 
(61). Unfortunately, even though the linearized operator L is simplified considerably, 
very little further analytical progress can be made and one has to turn to numerical 
computations. 

In order to compute the frequency co{k), for a given k, we need to calculate the 
linearized operator L. The computation proceeds as follows. First, we identify a basic 
vector Si which is a unit vector in six- dimensional space, with 1 at i-th component and 
otherwise. Then, we compute the matrix 

M{k) = (qi, . . . , qe) , q, = ( sin^ '^)JS, + J] L(s')^i • (78) 

s' ,m,k 

The frequencies u are then computed as generalized eigenvalues of 

M{k)S = uj\k)IoS. (79) 

In [Appendix C[ we consider a simplified pedagogical case when the polymer is only 
allowed to twist, which leads to (almost) analytic expressions for the linear stability. 
Unfortunately, in our case, the entries of 6x6 matrix M{k) for every wavenumber k 
have to be computed numerically as outlined above. It is not a difficult or numerically 
challenging computation though, and the computation of dispersion relation for several 
hundred values of k only takes a few seconds in Matlab on a standard desktop. The 
results of these computations are presented below. 
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8.2. Results of linear stability computations 

For the rod we are considering here, in the absence of the nonlocal interactions the 
rods are neutrally stable, as the elastic tensor J is diagonal with positive entries. 
It is therefore interesting that nonlocal terms introduce instability, corresponding to 



Im(ci;(fc)) > 0. Because u{k) only enters as uP in (79), the instability occurs when 
A = uP'{k) becomes negative. Physically, this instability is connected to the inclination of 
the rod to minimize its energy and properly align the dipole moments of each bouquet by 
twisting, as we have seen in the minimum energy calculation in Sec. |4j Mathematically, 
the instability corresponds to the eigenvalues of the linearization matrix M becoming 
negative for some k. However, one needs to keep in mind the discrete nature of our rods, 
making only certain values of k possible. Since we rescale the length by Iq, the distance 
between the centers of the bouquets, it is natural to take k G [0, 27r]. In the continuum 
case, there is no restriction on the wavenumber, so the condition for the instability is 
simply 

Instability = mino<fc<27r eigenvalues (M) < 0. 

On the other hand, for a discrete chain of length A^, k takes the values 2'Kn/N , where 
< 77. < A^. Thus, for a discrete chain 

Instability = minfc=27ri/Af eigenvalues (M) < 0. 

The difference between the discrete and continuous case is illustrated in Figure [7| The 
vertical axis shows the eigenvalue A = u'^{k) with A < corresponding to the instabihty. 
Only the part < < vr has physical meaning, since dispersion curves are symmetric 
with respect to reflection about k = n. There is an instability region for small k < k^,. 
However, a polymer of length = 9 (circles) does not have any allowed wavenumbers 
in the unstable region, whereas a polymer of length = 10 (blue squares) has one 
wavenumber /c = vr/S in the unstable region. Thus, from a physical point of view, even 
for a formally unstable situation, a short enough chain will be stable. There is some 
indication that such behavior is indeed observed in VDF polymers 



When the ionic strength / is increased, the Debye screening length decreases 



according to (26), thereby decreasing the electrostatic interaction. It is thus natural 
to assume that the rod's instability decreases for large values of I until finally the 
stabilizing elastic forces overcome the nonlocal forces. This is the case here. On the 
left side of Figure |8} we show the maximum unstable wavenumber k^ax as a function 
of the ionic strength /, and on the right side of this figure, we show the corresponding 
maximum stable rod length. We see that stabilization happens at / > 8.9M/1. 



9. Conclusion 



This paper has investigated the particular example of a simple molecule whose helical 
configurations possess the complex and intriguing energy landscape shown in Fig. |3] 
and Fig. |4j Yellow dots in these figures correspond to energy minima and thus show 
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Figure 8. The minimum unstable wavenumber (left) and its inverse, the maximum 
unstable wavelength, (right) are shown as functions of ionic strength in M/1. 



stationary conformations. We have derived a general scheme for analyzing linear 
stability of these states, particularly to elucidate the effects of torque on the molecular 
rod generated by non-local interactions of off-axis charge conformations. The stability 
analysis was facilitated by the 5* £^(3) symmetry of helical stationary configurations of 
the rod and it showed that non-local charge-interaction effects could induce instability of 
helical configurations due to the torques exerted on the rod by off-axis charges. This was 
illustrated by the instability of a linear polymer in its natural state. We have also shown 
how the increase of the ionic strength of the solution, in weakening the electrostatic 
interactions, leads to stabilization of the rod. The studies of general stability of stressed 
helical states will be forthcoming in our future work. 
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Appendix A. Definition of geometric properties of SE{3) group 

This Appendix describes the geometric structure of SE{3) group and defines the aspects 
of its adjoint and co-adjoint actions that are needed for the computations in the text. 
For more details and the theoretical framework, see [23]. 

Suppose G is a Lie group, and g and h are elements of G. Then, the AD operator 
- the conjugacy class of - is defined as 

ABgh = ghg-\ (A.l) 

for all g E G. Assume that h{t) changes smoothly with respect to a parameter t starting 
at the unit element of the group. Then, /i(0) = e and /i'(0) = r] is the velocity at the 
initial point, taken to be unity. Note that t] is an element of the tangent space at unity 
which is the Lie algebra of G. We denote this fact as T^G ^ q. In this notation, the 
Adjoint operation is defined as 



= ■ (A.2) 

t=o 

Note that Ad takes an element in rj E g and produces another element in q. Now, suppose 
g{e) is also varying with respect to a parameter e, and again g{0) = e, g'{0) = ^ G 0. In 
this notation, the adjoint operator ad is defined as 



ad^V = ^Adg(e)?7 



= ^V-V^=-- [^,V], (A.3) 

e=0 



where [■ , ■] is the commutator in the Lie algebra q. 

In order to derive equations, it is important to consider the co-adjoint operators 
Ad* and ad*. The operation Ad* : G x g* — )■ g*, is defined for g E G and a G 0* as 

{r] , Ad^a) = {Adgf] , a) (A.4) 

for every rj E Q, in terms of a suitable pairing (■,■): g* x g M. Similarly, the 
operation ad* : g x g* — )• g*, is defined for rj,^ E g and a E g* as 

{rj , ad^a) = (ad^r^ , a) . (A. 5) 

Let us now see how these formulas are expressed for the SE{3) group. The Lie 
group multiplication of two elements (Ai,ri) E SE{3) and {A2,r2) E SE{3), where 
Ai,A2 E 5*0(3) and ri,r2 E M^, is defined as follows: 

(Ai, r-i) ■ (As, ra) = (A1A2, Aira + r,) (A.6) 



Helical states of nonlocally interacting molecules 



29 



with the meaning of subsequent apphcation of rotation and shift. The inverse element 
is then 

(A,r)-^ = (A-\-A-V). (A.7) 

The tangent space T(Ao,ro)'S'£'(3) at a point (Aq, Tq) is defined as the space of derivatives 
of curves (A(t),r(t)) at t = given that A(0) = Aq. t(0) = Vq. In order to obtain the 
element of the tangent space at the identity - that is, the Lie algebra TcSE{3) ^ se(3) 
- we compute the derivative at Aq = Idsxs (a 3 x 3 unity matrix), Tq = 0. Hence, an 
element of this Lie algebra can be written as 

= (cj, v) = (uj, v) esc{3), 

where u denotes a skew-symmetric 3x3 matrix and v is a vector in three dimensions. 
Here, we may use the so-called "hat map" correspondence between the skew-symmetric 
matrices and vectors to define a vector a; G as Uij = —€^^00^, such that Cor = uj x r 
for all r. Thus, se(3) is a six-dimensional vector space, with the first three components 
having the physical meaning of the angular velocity, and the last three components 
being the linear velocity. 

Using this preliminary information, we are ready to define the adjoint actions. After 
some relatively straightforward computations, we have: 

AD(A,^)(A , r) = (AAA~\ -AAA'W + Ar + r). (A.8) 
Then if {u;,v)^^^{A{t),r{t))\t=o, 

Ad(A,r)(a), v) = (Aa;A~\ -Aa;A~ V Av) , (A. 9) 
and using ujA~^r — (jJ x A~^r — A~^ (Au? x r), 

Ad(A,r)(a) , v) = (AcuA"^ , -Au; x r + Av) . (A. 10) 

To express this in vector form, we note that for arbitrary u, Aa;A~^it = A(u; x A~^u) = 
Au; X It = (Au?) u, one has 

Ad(A,r)(^ , v) = (Auj , -Au; xr + Av) . (A. 11) 

Letting {u}i,cx.i) = ^(A(e), r(e))|e=o, one finds 

ad(^i,„i)(c^2 , 0:2) = i<^ii^2 - (^2(^1, -^20:1 + u:iOL2) , (A. 12) 

then because [a)i,C(;2]w = (u;i x uj^) x u for all tt, we can express this in vector form as 

ad(a,i,ai)(w2 , 0:2) = (wi X UJ2 , uJi x a2 - u;2 X ai) . (A.13) 

The physically relevant pairing between two elements (u?, a) e 5e(3) and (u, a) e 5e(3)* 
is given by 

((w, a) , (tt, a)) = u? • u a • a. (A.14) 

With this choice of pairing, we may also write 5e(3) as a pair of two 3D vectors. In that 
notation, the co-Adjoint operator is 

Ad(A,r.)-i (w, a) = (Au r X Aa , Aa) , (A. 15) 

and the co- adjoint action is given by 

ad*j^ ,^)(u, a) = (uxu; — axa, — u? x a) . (A. 16) 
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Figure Bl. (Supplementary, animated GIF file) Energy landscape for changing ionic 
strengths from / ~ IQ^^ to / = 10 M/1. All notations, axes and color scale are the 
same as in Figure [3j 

Appendix C. Twist dynamics of a straight polymer 

Appendix C.l. Simple case: twist dynamics of a linear rigid rod 

In this section, we consider the linear polymer drawn in Fig. [6] with the restriction that 
the charge bouquets can only twist about the axis, and only in the plane perpendicular 
to the axis. The rod itself is assumed to be completely rigid. This problem - in a 
slightly different configuration - was considered in [36] as a model of DNA dynamics. 
For convenience, we take the mass of the charges to be mo/2 so the moment of inertia 
has the value mo/o (no factor of 2), to coincide with the formulas derived in Sec. [s} We 
show two ways to analyze the linear stability of this problem. 

Standard solution 

The configuration space for this problem is described by the angles of rotation 0^. The 
coordinate for the positive and negative charges are given by 

rk,± = {kl, ±lo cos 4>k,±lo sin (pk)- 

The state 0, = for all i is an equilibrium state of the system. In order to linearize 
around that state, we proceed as follows. 

The distance between a charge at m-th unit and a charge at n-th unit depends on 
whether the charges are at the same or opposite sides of the chain. For the same side 
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we get 



dl 



ll {m - nf + P (cos (j)m - cos 0n)^ + P (sin 0^ - sin 



+ 



For charges on the opposite sides of the chain, 



■■\Jll{m- nf + /2 (cOS0rn + COS0^ 



+ P (sinf/)m + sine 



'/2 (m - n)^ + /2 



2T'n 2 



+ 



The electrostatic energy is positive for the charges on the same side, and negative for 
the charges on the opposite side, whereas Lennard- Jones energy only depends on the 
distance between the charges. Thus, the total potential energy P is given by the sum 



P 



\^J4\n^\^m,n,±U{d'^^) + U{d^^). 



The linearized equation of state is given by 



(C.l) 



(C.2) 



dt'^ d(j)m ' 

where / = uiqIq is the moment of inertia, mo being the mass of the charged particle. 
We choose the length to be measured in terms of /, and wavenumber < A; < 27r will 
be dimensionless. Assuming 0^ = e*'^'"~*'^*S', after some fairly simple algebra we obtain 
the following dispersion relation: 



/c.^(A:) = Jsin^^- 



U'{\m\) [/'(Vm2 + 4) 
\m\ + 4 



:i-e 



ikm\ 



m=— oo,m^0 

On summing up terms with opposite signs of m, one finds the real expression 

klo 



Iuj\k) = J sin' 



E 



m=l 



U'{\m\) ^ U'{^m? + A) 
\m\ ^Jm? + 4 



4sin^ 



mk 



(C.3) 



For the chosen values of physical parameters, a;' > and the twist dynamics is stable. 



Geometric method 

Assume that the axis of the rod is along x-axis, and all the charges in the undisturbed 
configuration are aligned along the z-axis. The deformations are rotations about the 
X-axis, which are given by the first coordinate in se(3) representation. The resulting 
rotation is also about the x-axis, so we need to look at the first row of the matrix M 



in (78). Thus, the stability of the problem of the twist about the axis is considered as 



simply the (1, 1) component of (79): 

Jnu;' = (M5i)i , (C.4) 

where we take 5*1 = (1, 0, 0, 0, 0, 0)"^ to be the vector of infinitesimal rotations about the 
X-axis. Again, after some rather straightforward algebra (not presented here) we see 
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that the right-hand side of (C.4) gives exactly (C.3). As verification of our code, we also 
compared the numerical results given by these two methods; they were identical within 
numerical accuracy of the computations. 

We note that (C.3) does not correspond to any eigenvalues, since the vector MSi 
is full, i.e. it has non-trivial components at other entries besides the first. These 
components arise because a realistic twist deformation about the x-axis induces twists 
and stretches in other directions. Nevertheless, we feel that such a simplified physical 
example still provides a useful verification procedure for the full stability calculation. 
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